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Abstract. We construct an auto-validated algorithm that calculates a close 
to identity change of variables which brings a general saddle point into a nor- 
mal form. The transformation is robust in the underlying vector field, and is 
analytic on a computable neighbourhood of the saddle point. The normal form 
is suitable for computations aimed at enclosing the flow close to the saddle, 
and the time it takes a trajectory to pass it. Several examples illustrate the 
usefulness of this method. 



1. Introduction 

It is well-known that computing a trajectory in the close vicinity of a fixed point 
is associated with many problems. Numerical integration schemes (silently) break 
down when the vector field tends to zero, and this usually results in completely 
inaccurate results. Indeed, as the norm of the vector field decreases, the fiow-time 
needed to pass a saddle increases without bound. This means that no integration 
scheme, rigorous or not, will function properly is this situation. There are, however, 
many instances where it is necessary to be able to follow the flow of a vector fleld 
arbitrarily close to a saddle. 

We present a completely automated, rigorous method that produces analytical 
estimates on the flow close to a given saddle. Equally important, it produces 
explicit bounds, for a given accuracy of the analytic estimates, on the size of the 
neighbourhood of the saddle on which the information is valid. This avoids the 
need to numerically integrate the flow near a saddle: once a trajectory comes close 
to the saddle, the bounds produced by our method give enclosures of where the 
trajectory exits the neighbourhood, and its associated flow-time. 

The approach is based on constructing a carefully chosen change of variables, 
that bring the original vector fleld into the robust normal form presented in [14lll5j. 
The present paper can be seen as a quantitative companion to [15], where several 
qualitative properties of robust normal forms are proved. Many of the ideas behind 
the algorithm can be found in [H], where they were used for estabHshing that the 
Lorenz equations support a strange attractor. In the present study we develop an 
algorithm for general planar real analytic vector flelds. 

Consider the planar vector fleld 

(1) i = Ax + F{x), 
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with A g <S, where S :— {diag(As, Au) : As < 0, Au > 0}, and where F is an analytic 
function, with F{x) = O(x^). Note that any vector field with a saddle fixed point 
can (locally) be brought into this form by an affine change of variables. 

The purpose of this paper is to describe, and implement, an algorithm that finds 
a square centred at the saddle in which we can enclose a trajectory and its fiow- 
time passing near the saddle. The output of the program includes estimates on the 
norms of the change of variables, its inverse, and the nonhnear part of the normal 
form, as well as the fiow-time for passing the saddle. 

2. Theoretical background and notation 

This paper addresses the algorithmic aspects of the planar case of the robust 
normal forms introduced in [14], and formalised in [15]. In order to simplify the 
formulae, we use vector and multiindex notation. The components of a vector are 
indexed by s and u for the stable and unstable direction, respectively. To make the 
presentation self-contained, we revise the necessary concepts from [15], but refer 
the reader to that paper for proofs and additional details. 

The structure of ([T]) implies that the stable and unstable manifold of the origin 
are tangent to the coordinate axes. Rather than attempting to find a coordinate 
change that completely Hnearises (HJ in accordance with Siegel's theorem [12], we 
compute normal forms that are robust in the sense that the set of eigenvalues where 
they exist is open and dense. This is crucial from a computational point of view, 
as we often only have an approximate knowledge of the eigenvalues. Our aim is to 
change ^ into the normal form 

(2) y^Ay + G{y), 

by an analytic change of coordinates, x — y + 4'{y)- We require that G, the non- 
linear part of the new vector field, is such that the invariant manifolds of the saddle 
are not only tangent to the coordinate axes, but actually coincide with them locally. 
We also require that the vector field is at least linear on these invariant manifolds. 
That is, if we let d{y) = min(|ys|, |2/ti|), then we ensure that Gi{y) = 0{d{yY), 
where I is the order of flatness. This means that if gm is a non-zero coefficient 
in the series expansion of G, then rag > I and ruu > I- We call the non-negative 
number |m| = niu + ms the order of m, and define the set W — {m G : \m\ > 2}. 
To formalise, let us split the space of multi-exponents into the sets 

V; {m e : rris < I 01 rUu < I}, 

U; := {m G : rng > I and m„ > I}. 

Now we can define the set of admissible linear parts of (JT)) that we consider: 

Ti := {A e S : m e V; mX - A^ ^ 0, i = w, s}. 

It is proved in [l5j that !Fi is open and has full Lebesgue measure in S. We will 
often use the notion of filters of a (formal) power-series: if f{x) = X]|m|>2 '-'^m-^™' 
we use the notation 

[f]ui = ^ ajnx"^, [fWi = ^ amx'^, and [/]„ ^ am- 

Also, we let denote the partial sum of the first d terms of /. We use the norms 
\y\ = maxdy^l, and ||/||r = max{|/(y)| : \y\ < r}. The r-disc is denoted by 
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Sr, and at times we use the notation A and A, to denote the eigenvalue with the 
smallest and largest absolute value, respectively. 

We are now ready to state the two main theorems from |15j : 

Theorem 2.1. Given an integer I > 2 and a system x = Kx + F(x) where F(x) = 
S|t?x|>2 o.mX^'^ is analytic, and A e Ti, there exists positive constants ro, ri, i^o, Ki 
and an analytic, close to identity change of variables x = y + 4>{y) with 

m\r<Kor^ (r<ro), 

such that X — Ax + F{x) is transformed into the normal form y = Ay + G{y) 
satisfying [G(y)](7, = G{y) and 

||G||, <Xir2' (r<ri). 

In the second theorem, we let ^' denote the flow oiij = Ay + G{x). 

Theorem 2.2. Under the same conditions as in Theorem \2.1\ and given any k > 
sufficiently small, there exists r > such that for any trajectory in 58^ starting from 
\xs\ = r, we have the following enclosure of its point of exit: 

'^u{y,Te{y)) = sign(?;„)r; 



r J \ r 

where Te [y) ( the exit time ) denotes the time spent inside of *Br • 

1 , r 1 r 
■ log -. — r < Te{y) < log ■ 



\u + K \yu\ \u-K \yu\ 

We will also use the following lemma from [15] . 

Lemma 2.3. If (As, Xu) are non-resonant for m G V;, then the divisors mX — Xi 
are bounded away from zero. Furthermore, for all orders \m\ > I + {I — I) 
we have the following sharp lower bound: 



\mX~Xi\ > |(|m| -0A+ - 1)A| 

Finally, the following lemma, which in principle appears in [14], will be used. 

Lemma 2.4. Ifr < ro(l — -ftToro), then (j) has a well-defined inverse, y — x-i-(j)~^{x) 
in \x\ < r* = r — \ \4>\\r, satisfying 

To prove the convergence of (j) and G we procede as in e.g. [fil [13], and use the 
method of majorants. If /, g : C" C", are two formal power series and \ fm\ < gm, 
for all multiindices m, and all the coefEcients of g are real and positive, we say that 
g majorises f, denoted by / -< 5. Thus, the convergence radius of / is at least as 
large as 5's. We will majorise in two steps; given some / : ^ C^, we construct 
g : — > C such that fi < g, for all z, and then construct /i : C ^ C such that 
g{z,z) -< h{z). 
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3. The algorithm 

In this section we will describe an algorithm that computes explicit bounds on 
the constants rg, ri, Kq, and n appearing in Theorems 12.11 and 12.21 This allows 
us to interrupt a numerical integration scheme of the original vector field JT]) , and 
instead use the analytical bounds from Theorem 12.21 to enclose the flow on passing 
the saddle, together with bounds on the time it takes to pass the saddle. 

The main ideas of the algorithm appear in [14], where robust normal forms are 
computed for the Lorenz system. In [14J, however, the algorithm was designed 
exclusively for that particular system. The purpose of this paper is to construct a 
general algorithm, that will take any planar vector field of the form |[T]), and trans- 
form it into |[2]), together with explicit bounds on the aforementioned constants. 

We note that several heuristic constants appear in the algorithm: t, rj, fi, p, 
Nq, ea, and e. In the actual implementation all of these can be set by the 
user in a configuration file. The algorithm has been implemented in a C++ program 
using the C-XSC package [51 [5] for interval arithmetic [U [U [H [TO] . For automatic 
differentiation [4] we use a modified version of the Taylor arithmetic package [2] . 

3.1. Outline. The algorithm has four main parts that will be described in detail 
below: 

1. we compute an I such that A e J^;, with I < l, where t is a user-provided 
order of flatness. 

2. we compute the flrst few terms of the formal power series solution of the 
functional equation for the change of coordinates x — y + 4>{y) using auto- 
matic differentiation. These first terms are used to estimate bounds on a 
majorant of (j). These estimates are then, finally, used to prove the analyt- 
icity of (j), using induction. 

3. we do the same kind of estimates for G; we compute some terms in the 
formal power series solution of a second functional equation, and use these 
to prove the analyticity of G by induction. 

4. using the estimates on the coefiicients of the analytic functions (j), and G, 
we estimate the constants Kq and k that enable the user to switch from a 
numerical integration scheme to the analytic estimates from Theorem 12.21 

3.2. Verifying A £ J^i. We want to determine I such that A e J^/. We will do this 
by first constructing V;, and then removing its members that cause resonances. 

Proposition 3.1. //, for i = I, I 

i N and i N, 

then A e Ti+i . 

Proof. For Z = 1, we note that Vi = {(0, i), {i, 0)}i>2. The potential resonances are 
given by 

'rriu^u — Ai = and m^As — Ai = 

and it is clear that no member of Vi satisfies any of these two equations. Hence, 
Ti =S. 

For / > 1, we have the recursive relation: 
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Thus we only have to consider the following potential resonances: 

'rriu^u + l^s — Ai = and ZA„ + m^As — Ai = 

with niu, rus > I. For the case z = s, we get 

rriu^u + (l ^ l)As — and IXu + {iris — l)Xs = 

with solutions m„ — {I — l)^^^ and = 1 + l^f-, respectively. Analogously, for 
the case i — u, we get 

(m„ — 1)A„ + IXs = and {I — 1)A„ + m^A^ = 

with solutions m„ = 1 + and nis = {I — 1)^^-, respectively. Therefore it 

suffices to enforce 

(3) i-^^N and i^^N, 

Au 

to estabHsh that A G Ti+i. □ 
It follows from Proposition 13.11 that we have the relation: 

Ti+i = Ti\{K e S : e N, or £ N}. 

To write a program that checks the condition in Proposition 13.11 is simple, and the 
algorithm returns a lower estimate on the largest I less than l, such that A e JT;. 

3.3. Computing and its radius of convergence. By inserting x = y + 4>{y) 
into (HJ), differentiating directly, and comparing the sides, we get: 

(/ + Dq^)y = A(y + 0(y)) + + 

By inserting into ([2]), and simplifying, we get: 

(4) D^iy)Ay - A0(y) = F(y + 4>iy)) - Dc^{y)G{y) - G{y). 
Let La be the operator 

La0 = D(t>{y)Ky - A(j){y), 

where we note that {Lj^{yY^))i — {raX — Xi)yY^ . 

Recall, we want to compute a normal form ^ which is /-fiat, that is [GJu, — G, 
and the non-fiat terms in which we want to cancel with 0, come from F. 
Therefore, by filtering on the component level, we get the following two functional 
equations for (pi and Gii 

(5) (iA0), = [F,(y-|-0(y))]v, 



(6) 



G, ^ my + <l>{y))]i 



^GM ~ —GM 



Since A e JT/, and [0]v; = (j) hy construction, we can solve ^ recursively. 

To bound the solutions of ^ we want to procede as in [14], and prove the 
convergence of the change of variables using majorants and induction. Two heuristic 
constants no, and ni > riQ are needed. They determine the range of coefficients of 
the formal power series of (p, that should be used in the induction proof. Let 



N{1) :=/ + 



(l-l) 



X 
X 
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be the constant from Lemma [2.31 from which the explicit lower bound holds. For 
the induction to work it is required that ni > N{1). 

We put ni = \{1 + ri)N{l)'] , and no = L^^(OJ , where 77 > 0, and -1 < ^ < 77 
are two given constants. 

Let 0i(a;) = J2'\^n\=2'^i,™,x"^ be the sought change of variables. We will com- 
pute the Qfi^m's with \m\ < ni using automatic differentiation, and then put ctk = 
S|m|=fe ™^^(l'^s,m|, |aM,m|)- The dfe's will be used as the first terms in a majorant 
of (j)s and Sometimes we will use di = 1, to simplify the argument of some 
functions. 

To calculate ai^m, with \m\ = k, we evaluate a fc- Taylor model of (a; + , 
and divide its mth term by mX — A, : 

I Am — Xi\ 

Note, the coefficients at a certain level only depend on the previous levels. This 
is because F does not contain constant or linear terms. 

If no and rti are sufficiently large, then the first terms computed above are a 
good approximation of a majorant 0, and we use these to determine an approximate 
radius of convergence for cj). The validity of this radius of convergence will be proved 
later. Therefore we determine, using a least squares estimator, constants C and M, 
such that 

dfe < cm'', no < k < ni. 

Thus, a candidate radius of convergence is s := -p, which needs to be verified. 
We will consider a slightly larger majorant of (pi. If 

00 

\m\=2 

we define 

Cfe := ^ max(|cs,m|, |c„,„|), 

\m\—k 

and set 

00 

k=2 

F is clearly a majorant of Fi. We define, 

(8) A:=± c.s'^-^ + ( "^-"" + "^""" (0 ' (P + 3) 

where p is a given natural number. 
Lemma 3.2. F{x) < A\x\'^ , on \x\ < s. 

Proof. The terms of F up to order p are clearly bounded by the left sum in 
since Cfc > 0, and \x\ < s. For the coefficients Ci^m, standard Cauchy-estimates give 
I < Thus, since there are (fc + 1) terms with \m\ = k, 
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Using C = 2s, this yields 

EZ.+iCk^" < m\\2s + \\Fu\Us)EZ,+Ak + i){^)' 



(9) 



- IP ^ l^k^p+lC^ + ^) [Ts 



I V fe-2 



+ Er=p+i(fc-i)G)'"' 



Let 

n{k) := (/c - Z)A+ (/ - 1)A 
be the lower bound on |mA — Ai| from Lemma [ 
Proposition 3.3. // 



□ 



then 4>i is analytic on Ss . 



k=l 



Proof. By the above lemma, (f)i is majorised by 0, where 0"^ = X]fcl2 ^kX^ is as 
above, and 

Note that no and ni are constructed so that max(2no, Ni) < ni. Thus, (for n > ni) 

A -A 
'^"+1 = 7v — r~rT 2^ afea„+i-fc. 

r2(n + 1) ^ 

Assume, for some n > rii, we have proved that dfe < CM'' holds for uq < k < n. 
If we can prove that ci„+i < CAf"+^, the convergence of (j) follows by induction. 

A f v^^o - - I \-^n—no 

""+1 = n(n+l) (^Z^fc^l Q^fcan+l-fc + Z^fc=no + l «fcQ!„+l-fc 

+ ELn-no + l afc^n+l-fc) 

- U{n+1) [■^ l^k=l aka„+l-k + 2^k=n„ + l "fcttn+l-fej 

(use the induction hypothesis, at < CM'', uq < k < n) 

< (2 ElU a^CM-^^-" + CHI-+^] 
= n(;^(2E^li«feM-'= + (n-2no)C)CM«+i 

< n^{2EZi^kM~''+nC)CM-+^ 

The expression before CM"+^ is decreasing in n for n > ni, since fl{n) ~ n|A|. 
Therefore, to prove the induction step, it sufHces to prove that the expression is 
less than one for n = ni. Thus, the close to identity change of variables converges 
to an analytic function on |x| < = s. □ 
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3.4. Computing G and its radius of convergence. It was proved above that 
there exists an analytic change of variables on \x\ < r^. After changing the 
coordinates, the system is of the form 

X = Kx + G{x), 

where [G]v, = 0. We want to estimate the radius of convergence of G. 

Let (j) be as in the above section, then a majorant for G is given by G, defined 

by 

(10) gk = [F{^ + ^''-') + 2{4>''+'-''yG''-%, 

where we note that gk = 0, for < fc < 2/, since [GJvj = 0. The g^'s are computed 
using automatic differentiation. 

We use ifTOj) to compute gk, for 21 < k < Nq, where Nq is an integer larger 
than 21. These values are used to compute candidate constants D and K, such that 
gk < DK^ . To do this, we again use our least squares estimator. We require that 
K > M; the reason being that will be used to estimate the radius of convergence 
for G, and G is only of interest within the radius of convergence of (j>. Let 

Hn) ■■= A[{2j:lUakM-' + C{n-2no))%{fr^') 

kauK^-' + CM (f ^""l^^y^ ' 

Proposition 3.4. //'I'(A^g) < 1, then G is analytic on Si^-i. 

Proof. Assume that we have proved gk < DK^ , for k < n. If we can prove that 
gn+i < DK"~^^, the convergence of G follows by induction. As in the proof of 
Proposition 13. 3^ we use the constant A to get a bound on F, which gives us the 
bound 

n n+2-2l 
(jn+l < A^^akO-n+l-k +2 ^ kakgn+2-k- 

k=l k=2 

We call the first sum Si, and the second sum S2- If we can prove that AT,i + 2E2 is 
bounded by 'i'{n)DK"^^ , where : N ^ R is a decreasing function, we are done. 



= (2 Efcli akM-'' + C{n - 2no)) CAf "+i 

< ((2 EZi e^kM-^ + C{n 2no)) % (#)"+^) DK-^\ 

since K > M, the bound on Ei is decreasing in n. 

^2 < D (j:Z2 kakK^+^-'' + C j:itl:+i kM'^K-+^-^) 

= (J:Z2 kakK'-' + CM (f (!^£±il^^^ DK-+^ 

the expression in front of DK^'^^ is independent of n. Thus, ^' is a decreasing 
function, which proves the bound gk < DK'' for all k. The analyticity of G on 
ri := i follows. □ 
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3.5. Computing the bounds. Let, tq = e^r^, r2 = ec^i, and = e min(ro, r2), 
where < e^, ec, e < 1, are given numbers. 

That is, is the domain of (j) that we will use to estimate ||(^|| and 
and *8r2 is the domain of G that we will use to estimate k. To ensure that our 
estimates hold we may never leave these domains. is the box where we will 
actually change coordinates. 

To guarantee that the change of variables is done in the domain of G, we need 
that + r^Ko < r2. By construction, the flow stays inside the domain of G, since 
the only place that the flow can leave the box is on the unstable side. To guarantee 
that the flnal change of coordinates is done in the part of the domain of </>, where 
the estimate on holds, we need that < ro(l — KqTo), since then r* < tq, 

where r* is such that = r* — \ \(t>\\r' ■ 

To compute Kq we note that on \x\ < rg, we have that 



Thus, if we put 



(11) i^o:=^d.rr+CM^(^''^°) 



then 



1 - 



To estimate \\(f> ^\\r, we need to find r* such that r = r* — ||0l|r* since then, by 
Lemma [231 < ll'/'llr*. A trivial calculation yield 



1 I 1 r 

— r — 



2Ko ■ Y4i^2 Ko' 
Thus, 



(12) \\cb-% < Ko{r*f < 



1 1 r 



— r — 



2Ko ' \i AKl Ko' 



Finally, the constant k is computed as 



We want that k <^ min(— A^, A„, |As + A„|); if this is not the case, we decrease r2 
and/or r^. 

4. Examples 

4.1. Example 1. We start with a simple example that also illustrates how the 
results depend on the distance from resonance. The vector field under study is 

(13) in = x^,+ 0.05x2 _ Q 95^^ 

+ ^((438.4905 - 25.2469a;s - 452.78992:2)a;„ - 741.0341a;3/3) 

which has previously been examined in [7]. It is a perturbation of a Hamiltonian 
system, given by <5 = 0. The Hamiltonian system has a resonance of fiat-order 
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1, since at a saddle of a planar Hamiltonian vector field the stable and unstable 
eigenvalues have the same modulus. 

We describe the results in detail for 5 = 10~^, and also include Table [H which 
illustrates how the convergence radii, and norm bounds depend on the distance 
from the resonance. We have chosen I = 10, since this is the lowest value of I that, 
after optimisation of — 0.1, ea = 0.5, and e = 0.9, yields k < 2~^^, which is the 
machine precision using IEEE double precision fioating point arithmetic. 

We start by introducing the linear change of variables, x = T^, that transforms 
(fT3| to the form |(T]) . Note that this transformation yields interval enclosures of the 
eigenvalues 

As = -0.77978852302649||, A„ = 1.21827902302649i, 

which are used during the computations. The diagonalised system is put into the 
algorithm. 

The algorithm starts by verifying that A g Tio, and computes 7V(10) — 25. 
Using fj, = 0.2, and t] = 0.08, yields no = 13, and ni = 30. Therefore, we need 
to internally represent all functions by their Taylor models of order 30. Next, 0^°, 
and (/)^° are computed using the recursive formula |[7|), and used to compute 
The least squares estimator of the coefficients of 0"^° yields 

dfe < 0.08 X 397'=, no < fc < ni, 

i.e. C = 0.08, and M — 397. We compute and use s — as the candidate 
radius of convergence to compute A — 1.64. To prove that (j) converges we verify 
Proposition [Q This yields = 2.52 x 10^^, vq = 2.52 x lO""*, and Equation [TT] 
gives Ko = 22.6. 

The algorithm now turns to the majorisation of G. We compute cjk, for 21 < k < 
4:1, using the recursive formula ifTOj) . The least squares estimator of the coefficients 
of G, yields 

gk < 1.03 x 10"^^ X 2490*=, 21 < k < 41, 

i.e. £1 = 1.03 X IQ-^^, and K = 2490*^. These values are used to verify Proposition 
[331 This yields n = 4.02 x 10"", rs = 2.01 x lO"'', = 1.81 x 10"", and 
K = 9.74 x 10-21. 

Finally, we verify that we compute within the domains of validity of the con- 
stants Kq, and K. When we enter the box |^| < ra, we apply (j), which al- 
ters the coefficient by at the most Kor^, that is we might start computing with 
Uu — {1 + Kor3)r3 < 1.82 x 10"''. Inequalitv 12.21 gives the bound, ys = ^{y, Te{y) < 

(1 + Kor3)r3 f (il^^iplllj ^" " < 1.84 x 10"^ We use EquationlH and compute 

< ITUs^V'^^ \J jk^ ~ — 8x10^''. Thus, the fiow exits the computations 

inside of the box 1^1 < 1.85 x lO"**, which is inside of |^| < min(ro, r2, ro(l — i^oro)) = 
r2 = 2.01 X lO^'', where the bounds on Kq, k, and are valid. 



4.2. Example 2. In our second example, we follow a solution curve close to a 
graphic. A graphic is an invariant set of a fiow consisting of saddles and separatrices, 
see e.g. Consider the following vector field 



X = {Sx + y){x^ - 1) 



(1^) y = {-x + 5y){y^-l) 
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K 




10- 


-3 


2.52 X 10- 


-3 


4.02 X 10" 


-4 


1.81 X 10" 


-4 


22.6 


9.74 X 10" 


-21 


10- 


-5 


6.85 X 10- 


-2 


1.35 X 10" 


-2 


6.09 X 10" 


-3 


6.52 


7.17 X 10- 


-19 


10- 


-7 


9.97 X 10- 


-3 


1.46 X 10- 


-3 


6.59 X 10" 


-4 


7.45 X 10-^1 


3.15 X 10- 


-20 


10- 


-9 


1.14 X 10- 


-3 


1.46 X 10- 


-4 


6.59 X 10" 


-5 


8.42 X 10-^2 


3.13 X 10- 


-21 


10- 


11 


1.30 X 10- 


-4 


1.46 X 10- 


-5 


6.59 X 10" 


-6 


9.67 X 10-^3 


3.14 X 10- 


-22 


10- 


13 


1.21 X 10- 


-5 


1.45 X 10- 


-6 


6.53 X 10" 


-7 


8.94 X 10+^ 


4.46 X 10- 


-23 


10- 


15 


1.07 X 10- 


-6 


1.28 X 10- 


-7 


5.79 X 10" 


-8 


7.90 X 10-^^ 


1.37 X 10- 


-22 


10- 


17 


8.62 X 10- 


-8 


1.11 X 10- 


-8 


5.02 X 10" 


-9 


8.28 X 10-^*^ 


3.08 X 10- 


-20 



Table 1. Convergence radii and norm estimates in Example 1 as 
S is varied. 




Figure 1. Phase portrait of the system from Example 2. 



where we will consider 5 = —0.2. If S ~ 0, this is a Hamiltonian field with the first 
integral H — ^ +x +y ^ There are five critical points, an unstable focus [centre 
if (5 = 0] at the origin and four saddles at (±1,±1), see Figure [H This example is 
simple enough so that we can determine most qualitative properties by hand, which 
allows us to focus our attention to the application of our algorithm. 

The curves a; = ±1 and y = ±1 are invariant under the fiow of In fact, 

they are the separatrices of the saddles. We only consider the fiow inside of the 
graphic. For S = 0, H = corresponds to the origin, and H = ^ to the graphic. To 
determine the properties of the fiow of lfT4|) . we begin by noting that, for 6 < 0, the 
vector field is transversal to the solution curves of the unperturbed system. Indeed, 
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let r2 = -^'"'+"'+^\ then 

2rf — —y^xx — yx^ij + yy + xx 

= {-y^x + x){Sx + y){x^ - 1) + {-yx^ + y){-x + 5y){y'^ - 1) 

= 5{x{x^ - + x)+ y{y'^ - l){-yx^ + y)) 

= S{x'^{-y'^x^ + x^) + y'^{-y'^x^ + Y^) - 2r'^ + x^y^) 

= (5(2r2(x2 + y2) _ - x^y^) ^ S{x^ + y^)i2r^ - 1) 

> 

Thus, if we leave the neighbourhood of one saddle on if = C, we enter the 
neighbourhood of the next one outside oi H = C . It follows that we do not need 
any numerical integrator to estimate the distance to the graphic from above. We 
will start in a neighbourhood of a saddle, use our computed analytical estimates to 
pass it and then enter at the next one at the same level curve of H. In addition 
we only need to consider one of the saddles, since the system is symmetric. We 
therefore translate (—1, —1) to the origin and get the system: 

^^^^ X — 2.2tZ/g X r^x g O.l^Z/g j 2x r^x g j 1.3.^^ 

Xii — 1.8x^ O.lx^ ~t~ x^Xs O.Tx^ 2xiiXs 

Our program yields the output shown in Table [2l 

Resonance of order 9 detected 
1=9, M_l = 19 

Order of Taylor approximations = 23 
n_0 = 10, n_l = 23 

C = [0.0180,0.0181] M= [3.9853,3.9854] 
A = [4.5520,4.5521] 

Phi is analytic on the disk with radius = [0.2509, 0.2510] 

K_0 <= [2.3223,2.3224] on the disk r_0 = [0.0376,0.0377] 

D = [1.2701E-010,1.2702E-010] K = [12.5616,12.5617] 

G is analytic on the disk with radius = [0.0796,0.0797] 

kappa <= [1.5544E-017,1.5545E-017] on the disk r_2 = [0.0262,0.0263] 

We recommend that you change to the normal form 

on the disk with radius r_3 = [0.0210,0.0211] 

Table 2. The output generated by the program in example 2. 



We will change to the normal form on So. 02, i-e. r — 0.02, and consider the 
trajectory that starts at (xs, a;„) = (0.02, 0.01) in the translated coordinate system. 
By Theorem [221 together with the bounds on Kq and k, we can calculate where it 
will leave So. 02- We start the calculation with < (1 + KQXi)xi, and then use our 
bounds on the flow inside So. 02, to get the following bound at y„ = 0.02 

2/s < (1 + Karjr I 

Thus, on the outgoing stable coordinate we have the following bound, 

(1 + KqXu)Xi, 



X. < (1 + ifor)r j +\\r% 

By the transversality and symmetry properties of the system, we enter the neigh- 
bourhood of the next saddle outside of (r, Xg), and so on. If we follow our trajectory 
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we get the upper bounds on its distance from the graphic, and lower bounds on the 
lap times, shown in Table [31 



lap 




laptime 





0.01 





1 


7.1 X 10-3 


1.7 


2 


3.0 X 10-3 


2.8 


3 


3.8 X lO-'' 


5.6 


4 


3.7 X 10-6 


12 


5 


1.2 X 10-1" 


26 


6 


3.0 X 10-1^ 


58 



Table 3. Converging to the graphic. 



To compare with the performance of a standard numerical integrator we do the 
same computations using the ode45 solver in MATLAB, which incorrectly starts 
fluctuating around Xu = 10-^, the result is shown in Table |4l 



lap 






laptime 





0.01 







1 


3.9 X 10- 


7 


23 


2 


8.0 X 10- 


8 


38 


3 


1.3 X 10- 


7 


39 


4 


7.9 X 10- 


8 


39 


5 


9.8 X 10" 


8 


39 


6 


1.3 X 10- 


7 


39 



Table 4. Numerical integration close to the graphic. 
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